"""
    Computational Economics
    3.2: More NumPy
    http://johnstachurski.net/lectures/more_numpy.html

    DEFINITIONS
    Linear Algebra
    Linear algebra is a branch of mathematics that studies vector spaces, also
    called linear spaces, along with linear functions that input one vector and
    output another.
    http://en.wikipedia.org/wiki/Linear_algebra

    Matrix
    In mathematics, a matrix (plural matrices, or less commonly matrixes) is a
    rectangular array of numbers, symbols, or expressions.
    http://en.wikipedia.org/wiki/Matrix_%28mathematics%29

    Determinant
    ...the determinant of a square matrix is a value computed from the elements
    of the matrix by certain, equivalent rules.
    http://en.wikipedia.org/wiki/Determinant#History

    Inverse
    If A is a square matrix and we can find another matrix of the same size,
    say B, such that then we call A invertible and we say that B is an inverse
    of the matrix A.
    http://tutorial.math.lamar.edu/Classes/LinAlg/InverseMatrices.aspx

    Matrix inversion plays a significant role in computer graphics,
    particularly in 3D graphics rendering and 3D simulations.
    http://en.wikipedia.org/wiki/Inverse_matrix

    Eigenvalue
    The eigenvectors of a square matrix are the non-zero vectors that, after
    being multiplied by the matrix, either remain proportional to the original
    vector (i.e., change only in magnitude, not in direction) or become zero.
    http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors#Applications

    Probability Mass Function
    In probability theory and statistics, a probability mass function (pmf) is
    a function that gives the probability that a discrete random variable is
    exactly equal to some value.
    http://en.wikipedia.org/wiki/Probability_mass_function


    REFERENCES
    http://en.wikipedia.org/wiki/Linear_algebra


    PROBLEM



"""


#
# EXERCISES
#
from random import uniform

def exercises():
    runs = 5

    q = [0.25, 0.75]
    for i in range(runs):
        print 'run %s: q=%s' % (i, q)
        print sample(q)


def test_prob_mass_function():
    """run 10000 times and compare results"""
    tries = 10000
    q = [0.1, 0.2, 0.3, 0.4]
    r = lambda f: '%.2f' % (f)

    counter = {}
    for n in range(tries):
        v = prob_mass_function(q)
        if v not in counter:
            counter[v] = 1
        else:
            counter[v] += 1

    output = []
    for n in counter:
        # (n, weight, expected, result)
        output.append( (n, r(q[n]), r(q[n] * tries), r(counter[n])) )

    print output


def prob_mass_function(l):
    """think weighted array"""
    a = 0.0
    U = uniform(0,1)
    N = len(l)
    for i in range(N):
        if a < U <= a + l[i]:
            return i
        a = a + l[i]


def sample(q):
    """
    Returns i with probability q[i], where q is an array
    (e.g., list or tuple).
    """
    a = 0.0
    U = uniform(0,1)
    N = len(q)
    for i in range(N):
        print "return i (%s) when U (%s) between a (%s) and a + q[i] (%s)" % (
            i, U, a, a + q[i])
        if a < U <= a + q[i]:
            print 'return i=%s' % (i)
            return i
        a = a + q[i]


#
# PROBLEMS
#
from numpy import searchsorted, cumsum

def problem_1():
    """
    Make the following two improvements to this code:

    1. Speed it up using NumPy, avoiding for and while loops
           * Hint: Use searchsorted() and cumsum()
    2. Implement as a class, where each instance stores its own value of q
           * Give the class a draw() method, which returns one draw from
             {0,...,len(q)-1} according to q
                 o Or even better, write the method so that draw(n) returns n
                   draws from q
    """
    q = [0.1, 0.2, 0.3, 0.4]
    p = ProbMassDraw(q)
    print p.draw(1000)

class ProbMassDraw:

    def __init__(self, q):
        self.q = q

    def sample(self):
        U = uniform(0,1)
        return searchsorted(cumsum(self.q), U)

    def draw(self, n):
        results = {}
        for i in range(n):
            r = self.sample()
            if r in results:
                results[r] += 1
            else:
                results[r] = 1
        return results


# Problem 2
import numpy as np
import matplotlib.pyplot as plt


def problem_2():
    """
    Problem:

    * Make the __call__ method more efficient using NumPy.
    * Add a method which plots the ECDF over [a, b]
          o Here a and b are method parameters
    """
    samples = np.random.uniform(0, 1, size=10)
    print 'sample:', samples
    F = ECDF(samples)
    print F(0.5)  # Returned 0.29
    F.observations = np.random.uniform(0, 1, size=1000)
    print F(0.5)  # Returned 0.479
    F.plot()


class ECDF:

    def __init__(self, observations):
        self.observations = observations

    def __call__(self, x):
        f = lambda n: n <= x and 1 or 0
        f = np.vectorize(f)
        return float(sum(f(self.observations))) / len(self.observations)

    def inefficient_call(self, x):
        counter = 0.0
        for obs in self.observations:
            if obs <= x:
                counter += 1
        return counter / len(self.observations)

    def plot(self, a=None, b=None):
        if a == None:
            # Set up a reasonable default
            a = self.observations.min() - self.observations.std()
        if b == None:
            # Set up a reasonable default
            b = self.observations.max() + self.observations.std()
        X = np.linspace(a, b, num=100)
        f = np.vectorize(self.__call__)
        plt.plot(X, f(X))
        plt.show()



#
# MAIN
#
if __name__ == "__main__":
    #exercises()
    #test_prob_mass_function()
    #problem_1()
    problem_2()
    print '%s: ok' % (__file__)
